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Abstract 

The stochastic LiouviUe equations are employed to investigate the combined signatures of chem- 
ical exchange (two-state-jump) and spectral diffusion (coupling to an overdamped Brownian oscil- 
lator) in the coherent response of an anharmonic vibration to three femtosecond infrared pulses. 
Simulations reproduce the main features recently observed in the OD stretch of phenol in benzene. 
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I. INTRODUCTION 



Recent experiments had demonstrated that two-dimensional infrared (2D IR) hneshapes 
can probe the picosecond dynamics of chemical exchange by observing coherence trans- 
fer in molecular vibrations through time- dependent spectral jump a^i^i^ . This spectroscopic 
technique^"^ is an optical analogue of 2D NMR commonly used to study slower (ms) chem- 
ical exchange processes^'^'SiiS. In this three-pulse experiment (Fig 1) the first pulse creates 
a vibrational coherence, whose evolution (free induction decay) during the first interval ti 
is related to the absorption lineshape by a Fourier transform. During the second interval 
^2 the vibrational frequency changes by complexation with the solvent. Finally vibrational 
coherence is again created by the third pulse and detected during the third interval ^3. The 
correlations of the lineshapes in the first and the third intervals provide information on 
chemical exchange taking place during the second interval. 

One experimentii looked at the OD stretching mode of phenol in benzene whose fre- 
quency changes between 2665cm~^ (free) and 2631cm~^ (complexed). The 2D lineshapes 
monitor the dynamics of complex formation and dissociatioi>i. For short delays (^2 ~ 100/s) 
two diagonal peaks broadened by spectral diffusion were observed at the linear-absorption 
frequencies. Off-diagonal cross-peaks are induced by chemical exchange for longer times 
(~ lOps). Similar effects were seen in the overtone region. The rate constants for complex- 
ation were directly obtained from the peak kinetics. 

Similar hydrogen-bonding dynamics were also reported for the CN stretch of acetonitrile 
in methanol^. The relative intensity of the two peaks attributed to free and hydrogen- 
bonded CN shows strong temperature dependence. The reaction coordinate profile in both 
ground and excited vibrational states was found to be very similar. The weak variation of 
the chemical exchange rates with changes in the vibrational state justifies the application of 
stochastic chemical exchange modelsii for vibration frequency fluctuations. 

In this paper we use the stochastic Liouville equations (SLE) to simulate the effects of 
chemical exchange^^i^^ and spectral diffusion in 2D lineshapes of an anharmonic vibration. 
This formalism was recently applied to model the signatures of conformational fluctuations 
in trialanine in water^^ using a Brownian oscillator coordinate model of fluctuations. The 
SLE were also used to compare the Brownian oscillator and the four-state-jump model of 
hydrogen bonding fluctuation on the photon echo of the OH stretch of HOD in DOD^^. 
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In section |n] we introduce the anharmonic vibrational Hamiltionian subjected to fluctua- 
tions due to a two-state chemical-exchange and a Brownian oscillator coordinate. In the two 
state jump (TSJ) model of Kubo and Anderson2iiLi2ii2i2£ bath-assisted (incoherent) chemi- 
cal exchange is described by kinetic equations. The model shows several distinct dynamical 
regimes depending on the motional narrowing parameter (ratio of the frequency shift and 
relaxation rate). Spectral diffusionji is incorporated by coupling to an overdamped Brow- 
nian oscillator coordinate^i. The 2D IR signals are defined in section IIIIl The stochastic 
Liouville equation for the TSJ model is introduced in section HVl and the 2D IR signals are 
calculated for the limiting cases of slow and fast fluctuations and for short and long-time 
delays between pulses. The combined effect of spectral diffusion and TSJ is presented and 
analyzed in section All lineshape regimes observed ini are qualitatively reproduced. 



II. THE MODEL 



We consider a single anharmonic vibrational mode described by the Hamiltonian: 

H = hnB^B + hjB^B^BB (1) 

where i?^ (B) are boson creation (annihilation) operators ([5,5"''] = 1). Both frequency 
(f2) and anharmonicity (A) are subjected to fluctuations described by 

A = Ao + Aid, + A2Q + Agfx.g. (2) 

Qq and Ao are the average values. Qi and Ai describe stochastic frequency modulation 
by chemical exchange represented by the TSJ mode^^^ (state up u {cz = 1) and down d 
{az = —1))- c^z is the Pauli spin matrix . 




Q is an overdamped Brownian oscillator (BO) coordinate causing spectral diffusion and 
coupled to the vibration through the parameters Q2, and A2. l^s, and A3 allow the spectral 
diffusion to have a different magnitude in the two spin states. We shall divide hamiltonian 
(Eq. (d))) into three parts H = Hsv + Hqv + Hqsv where Hsv includes the first two (spin- 
vibration) terms in Eq.(j21) which do not depends on Q, and Hqv, Hqsv include the third 
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and fourth terms in Eq.Q respectively. The dipole interaction with the electric field E{t) 

Hi^t = -fiE{t) (5t + B) (3) 

Three vibrational levels (Fig 1) are accessible in a three-pulse experiment: the ground state 
\g), the first excited |e) = B^g) and the doubly excited |/) = {l/\/2)B''B^g) state, with 
energies 0, Q and i7 + A respectively. The corresponding 9-component Liouville space basis 
set for the density matrix is denoted li^iy')) = Uju' = e,f,g. The dipole moment 

matrix elements are /ieg = /^e/ = v^yU- 



III. THE THIRD ORDER RESPONSE AND 2D SIGNALS 

We consider an impulsive four wave mixing process whereby the vibrational mode (Eq.(|H)) 
interacts with three impulsive optical pulses with intervals ti,t2, and calculate the response 
at ts (see Fig 1). Each interaction with the field can create or annihilate one vibrational 
quantum at a time (Eq®). The equilibrium distribution is 

|P(0)))« = Ip(0))5|0)q|^^)) (4) 

where \p{0))s is equilibrium spin state and |0)q equilibrium Brownian coordinate density 
given by Eq. ()D1|) . The space 7i is a direct product of the vibrational, spin and BO space. 

The polarization generated in this experiment is described by the third order response 
function^ 

s^'\h,h,h) = (^^J eih)e{t2)e{^^^ (5) 

where /i*-"-*^ = fi[B + B'^ , ^]. Summing over final states is represented by {{I\h = ms{Q\QTr 
where {I\s = (1,1), (0|q = / dQ and tracing is over the vibrational degrees of freedom 
Tr = {{gg\ + ((ee| + ((//|. The time evolution between pulses is described by the Green's 
function Q{t) of the SLE defined in the space Ti. 

The Liouville space state cannot change between pulses, so that Green functions are 
block-diagonal. The response functions can be separated into six contributions from different 
Liouville space pathways. 

Ri{h,t2,ti) = {I\Geg,eg{t3)Qee,ee{t2)Geg,eg{tl)\p{0))QS 
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Rii{t3,t2,tl) = {I\Geg,eg{t3)Qee,ee{t2)Gge,ge{tl)\p{0))QS 

-Rm(^3) ^2; ^l) = mGeg,eg {t3)Qgg,gg{t2)Qeg,eg{tl)\piO))Qs 

Rivih-ih.tl) = {I\Qeg,eg{h)Qgg,gg{t2)Qge,ge{il)\p{^))QS 

Rv{h-,t2,ti) = —{I\GfeJe{t3)Gee,ee{t2)Geg,eg{tl)\p{0))QS 

Rvi{t3,t2,ti) = -{I\GfeJeit3)gee,eeih)Gge,geitl)\p{0))QS (6) 

These may be graphically represented by the Feynman diagrams^^ given in Fig 1. For the 
stochastic models of frequency fluctuations considered here we have Ri = Rm and Ru = Riy, 
so that we only need to consider four independent pathways^^. The Green's function matrix 
elements Gfe,fe{t3) etc. are now matrices in the joint spin and BO space. 
We shall display the signal using the mixed time-frequency representation 

/•oo /"OO 

Ra,{u3,t2,UJi) = / exp{iUiti+iU3t3)Ra,{t3,t2,ti)dt3dti (7) 

Jo Jo 

Eq. (jH)) then reads 

Ri{uJ3,t2,UJl) = {I\Geg,eg{i^3)Gee,ee{'t2)Geg,eg{i^l)\p{0))QS] 

and similarly for other pathways , where 

POO 

Giuj) = / Git)e''^'dt 
Jo 

We consider the coherent nonlinear signals generated in the ki = — ki + k2 + ka and 
kii = ki — k2 + ks phase-matching directions. Invoking the rotating wave approximation 
the ki (photon echo) signal is^ 

Sl{uJ3,t2,UJi) = ^0 {p'^g[Rii{u3,t2,l^l) + Rivi(^3,t2,UJl)] + plgPlfRvi{uJ3,t2,UJi)] (8) 

and for the kn signal we have 

Sll{uJ3,t2,UJi) = {p'^^g[Ri{u3,t2,iJl) + Riii{uJ3,t2,UJl)]+ plgPlfRv{uj3,t2,UJl)} (9) 

We shall also display the following combination which shows absorptive peaks^^^. 

Sa{uJ3, h, t^i) = -Im [^/(cjs, ta, -^i) + ^'//(cjg, ts, ^i)] (10) 

When during the interval ^3 the system has no memory about its state in ti the third 
order response functions are factorized into products of the linear response functions K{t) 

R^{t3MM) = K{t3)K{ti) 



Riiih, t2, h) = K{h)K*(ti); Ru{u3, t2, uji) = K{uj^)K\-uj{)- (11) 

Neglecting the overtone contributions and Rm it follows that the Sa signal is given by 
the product of the linear lineshapes Wa 

HSaIoos, uj,) = WaMWaM. (12) 

IV. STOCHASTIC LIOUVILLE EQUATIONS FOR THE TWO-STATE-JUMP 
MODEL 



We first consider the TSJ model^*i^^ by neglecting the spectral diffusion (setting 
A2,A3,fi2,^3 = 0). The total density matrix p has 18 components |z/z/'s)) given by the 
direct product of 9 Liouville space states and the two spin states, s = u,d. The 

SLEii-2& then reads 

dp 



dt 



Cpit). 



(13) 



The Liouville operator C is diagonal in the vibrational Liouville space variables, and is thus 
given by nine 2x2 diagonal blocks in spin space. 



L J vv's,v\v'^s' 



ss' 



c 



sv 



J VV' SMV'S 



(14) 



(15) 



where Ls describes the two state jump kinetics and Csv represents the coherent vibrational 
evolution that depends parametricaly on the spin state. 

Chemical exchange is described by the rate matrix Cs for the spin: 

kd ku 

The up (down) jump rates {kd) are connected by the detailed-balance relation ku/kd 
exp /9(erf — e„), where — e„ is the energy difference between the u and d states. 
The equilibrium density matrix is 

|p(0))) = \9g))\pms; |p(0))5 = TT^iT 1 ^" 

l^u + l^d \ kd 

We next turn to the vibrational part Csv = —{i/fi)[Hsv, ■ ■ ■]■ For the \gg)), |ee)), |//)) 



(16) 



blocks 



Csv 




Csv 




Csv 




99,99 




ee,ee 





c 



sv 



e9,eg 



= 0. The other blocks of Csv are: 

//,// 

-i{Qo + Qi) 

-i{no-^i] 



sv 



faJa 



-i(no + Ao) + Ai) 

-i{no + Ao) + + Ai] 

-^(2^0 + Ao) - i(21]i + Ai) 

-i(21]o + Ao) + z(2f]i + Ai; 



C 



sv 



C 



sv 



J v'v,v'v 



The remaining blocks are obtained by taking complex conjugates 

These matrices may be represented in a compact form by defining the energy heo and 
splitting hei parameters ej,^^ = 0, ef^ = fio, t\j^ = 2Qo + Aq, e[^^ = 0, e'f^ = Qi, e[^^ = 
2^1 + Ai. We then have 



sv 



-lie 



^0 y 



ne 



1 ^1 y 











-ne, 



^0 , 



(17) 

The Green's function solution of Eq. ()13|) is given by 2 x 2 blocks for each vibrational 
state of the density matrix 



exp (rjit) + 



exp (772^) 



^V2-Vi V2-V1 J \V1-V2 V1-V2 

where rjj are the eigenvalues of each block of C 



(18) 



Vi = 7; Heo - ^0 



en ' ) + A/ ^-^^±^ - (e^) - ei^y + .(eS^) - er^)(^. - K) 
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For the gg,ee and // space Csv = 0, and Eq.lfTHj) reads 

QeeAt) = Q^aAt) = Offjfit) ^ exp [Cst) = i + 

'(19) 

The linear response is given in Appendix 1X1 Closed form expressions for the various pathways 
are given in the Appendix El Below we discuss the 2D signals for limiting cases. 

In the slow bath limit Qi >> k (In Figs 2-6 k = = ky) no spin jumps occur during 
the ti,and ^3 intervals. In Fig 2a we display the 2D lineshapes Sa (Eq.^^). In all figures 



we give frequencies with respect to Qq (i-e- we set Qq = 0) and normalize the signal to have 
maximum absolute value of 1. For fct << 1 no spin jumps occur during the t2 interval and 
we only see 4 peaks (Fig 2a) two diagonal at Ui = = Qq ± Qi and and two overtone at 
{ui, UJ3) = {Qq + Qi, Qq + Aq + Qi + Ai) and {Qq — Qi, Qq + Aq — Qi — Ai) In the opposite 
kt2 >> 1 limit we see 8 peaks (Fig 2b). In addition to the previous four peaks cross peaks 
appears at {lui, lus) = {^0 + ^0, ^0 - Aq), {^0 - Aq, ^Iq + ^0), (^0 - ^^i, ^^o + + fii + 
Aq), (f2o + ^1, — ^1 + ^0 — ^i) due to the spin jumps. Diagrams i, ii, Hi, and iv describe 
the g to e peaks at 001,003 ~ Q. Diagrams v and vi show the e to f peaks at cji ~ + A. 

Note that the relation (Eq. ()B5jl ). similar to Eq. lfT^ is valid for individual peaks in the 
short time regime since peak shapes are not connected with any inhomogenity which carries 
memory, however Eq. (jl2j) does not apply for the whole 2D shape, since the memory of spin 
state erases the cross peaks. In contrast, for long ^2 (Fig 2b) this memory is lost and Eq. 
flT^ applies (excluding the overtone contributions). 

In Fig 2c we plot Sa (EqlTO|l in the fast ku,kd » fli limit. We see a fundamental peak 
at 001 = 003 = Qq and an overtone peak at ooi = Qi, 003 = Qo + Aq. Linear response in fast 
regime shows one peak (fiuctuational narrowing^''). Any memory of ti interval is lost in the 
^3 interval and the factorization fEqlTT|) is valid resulting in Eq. which is observed in Fig 
2c. 

V. TWO-STATE-JUMP WITH SPECTRAL DIFFUSION 

Spectral diffusion is described by a dimensionless overdamped Brownian oscillator coordi- 
nate Q (Eq.(^) whose variance 1 and its dynamics is given by the Fokker-Planck operator— 

where A is the relaxation rate. For fast fluctuations ^2/-^ << 1, spectral diffusion may 
be accounted for by simply adding dephasing rates fi^/A to the lineshapes. The complete 
solution of the Fokker-Planck equation is required for arbitrary fluctuation timescales. 

We define e^^^ = 0, e^'^ = ^2, = 2^2 + A2, and e^^^ = 0, e^'^ = ^3, e^^^ = 2^3 + A3 
The coupling of Q to the vibration is given by the Liouville operator Cqv 

^SQV 



- ^(^l:^-er)QSuuA'.iWzU (21) 



8 



The complete Liouville superoperator describing both TSJ and BO fluctuations is flnally 
given by: 

^tot = + ^SV + ^QV + ^SQV + (22) 

where Cs is given by Eq. (|THjl . Csv by Eqs. (|T7|l . Cqv^C-qsv by Eqs. fpTjl . and Cq by Eq. pn|l . 

The response function for the spectral diffusion model alone^^ (neglecting the spin, setting 
VLi = 0,Ai = 0) can be calculated using the second order cumulant expansion^^ and is given 
in Appendix O 

We next turn to the third order response (Eq.®). The equilibrium distribution is 

Ip(0))q5=|p(0))5|0)q (23) 

where [^(O))^ is given by Eq. (fTH|) and |0)q by Eq. ()D1|) and the flnal averaging should run 
over all degrees of freedom. 

{I\qs={I\s{0\q 

The Green's function matrix elements Gfeje{t-i) are now matrices in the joint spin and BO 
space 



iuj + C 



tot 



-1 



(24) 



Q{uJi), QiuJs) are calculated in Appendix [E] by expanding them in the Fokker-Planck 
eigenmo deai^^iSiSLSa. 

During ^2 the evolution is in the ee, gg, f f space where Cqs = and the Green's function 
may be factorized as 

Gee,ee{t2) = Qgg,gg{t2) = exp (^Csh^ exp (£Qt2) (25) 

where the flrst factor is given by Eq. ()19j) and the second term is the propagator of Fokker- 
Planck equation^i expanded in its eigenmodes fEq. ()Dl|) ) 



[exp (Cqt)]^!^ = Saf3 exp (-aAt) (26) 

For slow chemical exchange compared to the spectral diffusion {k^jka << cr, A) the TSJ 
peaks are well resolved and their lineshapes are determined by the spectral diffusion. Chem- 
ical exchange occurs on a much longer timescale and affects the cross peaks. 

We start by setting ^3, A3 = and examine fast and the slow spectral diffusion. The fast 
(motional narrowing) limit (A >> ^2) gives a Lorentzian absorption lineshape. Since the 



bath has no memory, peak shapes are less sensitive to ^2 and cross peaks due to chemical 
exchange only appear as ^2 is increased. Effects of fast fluctuations may be accounted for 
by adding dephasing rates 2D-lineshapes Sa (Eq- (llUp shown in Fig 3 are similar 

to the TSJ model (Fig 2), both are given by a product of Lorentzian lineshapes (in 1^1,^73). 
The important difference is that when accounting for fast spectral diffusion, the linewidth 
is determined by both exchange and dephasing rates {kd + h)~^^ i.e. it is different from 
the cross-peak timescale. In contrast the linewidth of TSJ model must be equal to exchange 
rate k obtained by analysis of the cross-peaks growth. 

Sa in the static A << f22 limit (Fig 4) shows Gaussian linear lineshapes. Static disorder is 
eliminated in the photon echo ki experiment. 2D peaks are elliptic with different "diagonal" 
and "off-diagonal" linewidths (Fig 4a)'^°. The intermediate regime may be observed at 
A^^ <t2 < k~^ when the bath loses memory (Fig 4b) and the 2D peaks become symmetric. 
For long times kcit2 >> 1 we see cross peaks induced by chemical exchange. Note the circular 
contours of Fig 4c, compared to the star-like contours for fast fluctuations (Fig 3). These 
may be understood by noting that the product of two Gaussians (which represent the linear 
lineshape for the slow case) is rotationally invariant, unlike the Lorentzian lineshapes in the 
fast limit. 

Fig 5 shows the ki and kn signals (Eqs. (jHl), and Q) • The ki contribution dominates 
at short times, and shows effective rephasing (photon echo) and the final 2DIR is similar to 
this dominating kj contribution (See Fig 4a,). The bath loses memory of its initial state in 
the third interval. For longer t2 » A^^ (right panels), the kj and kn contributions become 
comparable and 5*^ becomes symmetric (see Fig 4). 

We next explore some more general cases. First we take A2 7^ 0. The overtone peaks 
may then show different widths along and uJi. (overtone photon echoes are expected at 
^2^1 = (f^2 + ^2)1^3 ) Note that we have a slow bath (A < ^2) for the eg and a fast bath for 
the ef transition (A > ^2 + A2). These lineshapes are shown in Fig. 6. 

The inclusion of ^3 does not require additional numerical effort. This parameter allows 
the u and d peaks have a different width due the spectral diffusion. 

To demonstrate these effects and relate them to recent experiments^^ we have simulated 
the spectral diffusion using a single overdamped Brownian coordinate allowing a different 
width for u and d peaks. We used the splitting 2Aq = SAcm'^ (i.e. ~ l.Olps"^) and the 
exchange rates fc„ = 0.1ps~^ kd = 0.125ps~^ which were determined ini from the cross peaks 
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growth. Three regimes similar to those of Fig 4 were observedi suggesting that the bath is 
not fast. The intermediate regime, where the memory of the BO coordinate is lost (circular 
lineshape) but the cross peaks are still weak, was found on the 2ps timescale. We thus 
assumed for the relaxation rate A ~ OAps~^. Knowing A, Qi and are connected to the 
peak linewidth in linear spectra and can be estimated using the Fade approximant of a 
2-level syste m^^i^-^ . We have simulated the absorption spectra in Fig 7 and found lineshapes 
similar to experiment for ^2 = 0.33ps~^, = — 0.07ps~^. This completely determines the 
model for 2D IR spectroscopy (neglecting dynamics in the overtone), we have no additional 
free parameter. In Fig 18 we show that the predicted 2DIR lineshapes of the same model are 
in qualitative agreement with experiment. We see all three regimes, rephasing eliptic shapes 
, the relaxed Brownian oscillator with circular shape and chemical exchange cross-peaks at 
the proper timescales. The peaks have also the correct relative intensity the lower frequency 
peak is weaker but broader. 

In conclusion, we have demonstrated that the SLE may be used to model chemical ex- 
change in 2DIR spectroscopy. Spectral diffusion can no longer be accounted for in terms 
of dephasing rates, when its timescale may be observed in experiment and its proper de- 
scription must be combined with the two state jump which describe chemical exchange. The 
high temperature overdamped Brownian oscillator model for spectral diffusion with arbitrary 
fluctuation timescale reproduces the most significant features of recent experiments^. 

In contrast with calculations based on cumulant expansions at finite temperature o^^i^^i^^ 
our high temperature bath does not respond to the state of system, and its evolution in 
the ground and excited state is same. This may limit the applicability of SLE model when 
signatures of the vibrational Stokes shift are observed in 2DIR signals, however it simplifies 
the calculations. The lack of phase factors during the ^2 interval reduces the dynamics to 
classical level and allows large scaleMD simulations of environmental dynamics. 
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APPENDIX A: THE ABSORPTION LINESHAPE 



The linear response function is given bj» 



S^'\t^) = 9it^)-{{I\^,(-^git^)^,(-^\pm)n 



We define the contribution 



i^(tl) = (/|^e,,e,(tl)|p(0)) 



QS- 



This response function is connected to absorption hneshape 



W{ui) = i4gh-^Re I exp {iLUiti)K{ti)dti 



For the TSJ model (1^2,3 = 0) the result of Kuboi^ is recovered: 



(Al) 



W{uj) 



{kd + ku) (a; - fio - ^^i)' + [{^ - ^o){kd + ku) + n^ikd - 
For the spectral diffusion model we set Qi = and the linear response function is given 
by the second order cumulant expression (see Eq. ()C6|) ). 



K{ti) = exp {-geeiti)) 

APPENDIX B: RESPONSE FUNCTIONS FOR TWO-STATE-JUMP SPIN DY- 
NAMICS 



In the frequency domain the Green's function (Eq. ITH|1 matrix elements are: 



iuj + jC 



-1 



- 4^^^ + 4^^'^)^ - {^['^ - + - + et\kd + k^) + ^(e^) - et\kd - k^) 



Kd Kd — ^[^ — + ^0 ~ ^1 + ^1 J 



(Bl) 



Ri and Ra (EqEI) are given by: 



LJi ^ Qq)'^ + =F ^o){kd + ku) ± i^iikd - ku) 



12 





(a;3 - Oo)^ -^1 + i{uJ3 - ^o){kd + K) + i^i{kd - K) K + ka 

X |(A;„ + kdf - [{cus - ftoXcui t ^q) ± ^1] (ku + ka) + [{cui ^ llo)f^i ± (^^3 - {ka - k^) 

-i [{u;3 -no + u;iT ^o) {ku + k^f - (l^i ± Oi) {kl - kl)] 

where the upper (lower) sign is for Ri {Ra). Since the evolution in the excited |ee)) and in 
the ground \gg)) state is the same for stochastic jumps (where the bath is not affected by 
the system) we have 

Riii{i^3,t2,UJi) — Ri{us,t2,UJi); Riy{uJ3,t2, LOi) — Rii{L03,t2, LOi) (B2) 

We further have for Ry and R^f 

Ra{uJ3,t2,i^i) = - j—^TFTv?. 7^ 



X 



{(Ji =F ^^o)^ - of + i{u;i ^ ilo){kd + ku) ± i^i{kd - ku) 
1 

(u;3 - Qo - Ao)2 - (Qi + Ai)2 + {(003 - Qi - Ai)(kd + ku) + ^(^ll + Ai){kd - ku) 



^{(/^u + kd)"^ - [{003 -Qo- Ao)(a;i T ^^o) ± (^^i + ^1)^1] {ku + M 
+ [{ui ^ Oi)(Oi + Ai) ± {u3 -flo- Ao)fli] {kd - ku) 



[{uj3 - fio " Ao + cc^i T ^o){ku + k^)^ - {n, + Ai ± n^){kl - kl)] 

±m{^^+A,)k.kd '-'''''lr['l^'-^'^h 

r^d ~T Ku ^ 



where the upper (lower) sign is for Ry {Ry 



1. High temperature spin dynamics 

In the high temperature limit we set ku — kd = k. We first analyze the Ra contribu- 
tions. Here 

a.„,(t) = e-"+*"^°''<"*'"'^"'^''^* I (B3) 

\ ^ sinh {-qt) cosh {r]t) + ^ sinh {iqt) 

and 

Qge,ge{t) = Qlg,eg{t) 

13 



where rj = ^/k"^ — When > k we use analytic continuation (e.g. cosh{ix) = cosx 
and sinh [ix) = i sinx) and after some algebra we get for a = i, iiM 

g-A;(ti+t3)g-if^(i3±ti) 
Raits, hytl) = 



P\ f k'^\ 2k 

1 + — cosh [r]{ti + ta)] + 1 - ^ cosh [r]{ti - ts)] + — sinh [r/(ti + tg)] 



1 



{uJs-fto?~fll + 2i{u3-Qo)k 

-2i{uj3 -rio + uJiT ^o)k ±nl[l- exp (-2H2)] | 

We further have Rm = Ri, Ri^ = Ra. For a = v,vi we get 

1 

Rai^3, h, 1^1) 



1 



{003 -^0- Ao)2 - {n, + Ai)2 + 2t{uj3 - - Ai)A; 
X |4A;2 - [(a;3 - fio - Ao)K ^ ^]o) ± (^^i + ^1)^1] 
-2% [(cu3 - - Ao + cji T ^]o)A;] ± + Ai) (1 - exp [-2^2]) } 

where the upper (lower) sign corresponds to R^ [Rvi)- 



^ g-fc(ti+2fe+t3)g-in(t3±ii) ^ [cogh b(ti + tg)] - cosh [riiti - tg)]] (B4) 

where the upper (lower) sign stands for Ri (Ra). Eq. ((7|then gives 

1 



2. Slow spin dynamics 

When assume Qi » ku, kd no spin jumps occur during the ti,and intervals. Then 



kd + ^iei''^-el:'^+e[^^-e[''^-co) 



k^ + - - e['^ + e[''^ - u) 

The Green's function for ^2 is still given by Eq. (jl9p . 

In the vicinity of the peaks (Fig 2a) at e.g. Ui ^ VLq + Vti^ oj^ ^ VLq — Vti we have 

1 1 

i?i(w3,0,u;i) oc 



ku - i{uJi - r^o - ^1) kd - i{uj3 - ^0 + ^1) 
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{ki + {uji -no- ^^i)2)(A;2 + (cus-no + n 

Rii{uJ3,0,UJi) oc 



ku - i{ooi + fio + ^i) kd - i{uj3 - + Qi) 
_ kdku - {uJi + Qo + - Qq + Qi) + i [kdjuJi + Qp + Qi) + k^juj^ - + ^^i)] 

{kl + {UO, + fio + ^lYWd + (^3 - ^^0 + ^lY) 

In the absorptive signal (Eq. ljlOj) ) the dispersive parts (the second term in nominator) cancel 
out and the signal is given by simple lorentzian linear response peaks 

3. Fast spin fluctuations; motional narrowing 

In the kui kd» Vti limit we get 

^ _ -1 Iku- i{uj - VLq) ku 

''''' - ^o)' + ^o) {kd + ku) y kd kd- i{uj - fio) 

Qge,ge{^) = Qeg,eg{ — ^)* 

The ine shapes is insensitive to t2 delay one can set 

The response functions may now be factorized into products of the linear response functions 
(Eq. (HID). 

APPENDIX C: THE CUMULANT EXPANSION FOR THE SPECTRAL DIFFU- 
SION MODEL 

Neglecting spin dynamics by setting Ai = A3 = 0, i7i = = 0, the response functions 
may be calculated by the second order cumulant expansion^SiS which represents stochastic 
Gaussian fluctuations. 

The response function may be expressed in terms of the following four-point correlation 
functions^: 

F2{n,rs,r2,n) = {B{n)B{rs)B^{T2)B^{n))/2. (CI) 
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We then have^"^^: 



Rii(t3,t2,ti) = Fi{0,ti + t2,ti + t2 + t3,ti) 

Riii(t3,t2,ti) = Fi(ti + t2 + tsjti + t2,ti,0) 

Riv(t3,t2,ti) = Fi{0,ti,ti + t2 + ts,ti + t2) 

RvihMM) = F2{h,h + t2 + h,h + t2,o) 

Rm{h, t2, h) = F2(0, ti+t2+ h, h + t2, h) (C2) 



Here: 



-^"1(^4, n, T2, Ti) = exp[i{-QorA + ^on - ^IoT2 + QqTi) - /i(r4, T3, T2, n)] (C3) 

with: 

/l(^l, T2, Ts, T4) = geeihl) + QeeiUi) + gee{h2) + gee{Ul) - 
~5'ee(^3l) ~5'ee(^42)- 

where Uj = Xj — Tj and 

^2(^-4, i"3, ^2, ri) = exp[i{-non - i^o + ^o)n + {^o + ^o)t2 + ^on) - f2{n, Ta, rff^A) 
with 

/2(ri, Ta, Ta, T4) = fi-eefel) + fi-Z/fe) + Qeeitis) - fi'e/(t2l) - 

-9ef{h2) + fi'e/(t3l) + 5'ee(t32) + QeeiUl) - 

-gee{hl) - 9eeiti2) " fi-Zefe) " 5'/e(i43) + 5'/e(^42)- (C5) 

(yfafe are the hne broadening functiona^^iSi 

gee{t) = [exp (-At) + At - 1] 

Qefit) = (^]2 + A2)fi2 [exp (-At) + At - 1] 

^//(t) = (1^2 + A2)' [exp (-At) + At - 1] (C6) 

Note that Ri = Rm, and Ru = Riy as expected for stochastic models of frequency fluctua- 
tions . 
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APPENDIX D: SPECTRUM OF THE FOKKER-PLANCK EQUATION 



The eigenvectors of the Fokker-Planck operator Eq. ()20|) with eigenvalue — aA are given 

----- 



where is the Hermite polynomial 

HJx) 



-ire"" 



-e 



The matrix representation of bath densities and evolution matrices refer to to the basis 
{la)^}. The matrix elements of Cq are 



(D2) 



We use the recurrence relation: 



The Q coordinate is then represented by the tridiagonal matrix (Eq. ()D4I 



(D3) 



[Q]l3a — l3V25p^a+l + y^^''^'""^ 



(D4) 



APPENDIX E: MATRIX CONTINUED-FRACTION SOLUTION OF THE 
FOKKER-PLANCK EQUATION 



The complete Liouville superoperator is a 18 x 18 matrix in the joint vibrational and spin 
space. 

The Q variable is tridiagonal in the Fokker-Planck eigenbasis and the complete Liouvillean 
may be thus arranged in the tridiagonal block structure in the Brownian coordinate variable 



/ So ... 

Qi Qi Qt 




(El) 
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For Eq.(I22I) the blocks are: 



Qn = ioJ + ts + jCsv - nAl 

Qn = ^^QS 

The SLE may be solved in the frequency- domain using a matrix continued fractioni^^^LSS 
The Green's function is given by the inverse fEq]^^ of the tridiagonal matrix Eq. ()El|) . 
Starting with 

iu + C G{uj) = —1 
For the off- diagonal n ^ m elements we have 

QnOn-l,m + QnGn,m + Q^Qn+l,m = (E2) 

The diagonal n = m elements are 

QnGn-l,n + QnGn,n + Qn'5n+l,n = "1 (E3) 



The recursion relation Eq. ()E2|) is independent of the index m. Consequently we can 
introduce the matrices S^^S^ 

Gn±l,b = <SnGn,m (E4) 

Using Eq. ()E2|) these matrices may be solved iteratively. 

= Q Ql±l (E5) 



Combined with Eq. ()E3|) we obtain for the diagonal terms: 

^"'^ ^ QnS;. +Qn + Qtst ^^^^ 

while the off-diagonal Q{uj)nm are obtained from Eq. ()E5|) and ()E4|1 . For n > m 

Gn,m{^) = Sn-li^)<Sn-2i^) " " " (E7) 

and for n < m 

Gn,m{^) = 5~+i(^)5;+2(w) ■ ■ ■ S-{uj)gm,m{^) (E8) 



The full Green's function g{uj) can be calculated by using Eq. ()E5|) to find the connection 
matrices S^. We note that Sq is zero and truncate the recurrence relation for at some 
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level n by setting iS^ equal to zero. The matrix elements Q{uj)mm. can then be obtained 
using Eq. ()E6|) . All other matrix elements Q{uj)nm are calculated from Eq. ()E4|) . 
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Figure captions 



Fig 1 Top: Pulse configuration for a tliree pulse cohrent experiment 

BOottom: Feynman diagrams for the third order coherent response of an anharmonic 
vibrations. The diagrams correspond respectively to the six terms in Eq. (jU)). 

Fig 2 (Color Onhne) 2D signals Sa {Eq.^ for the TSJ model (Eqs. (HI, (UHl)), (a), static 
limit ^i/ku = 5 for short delay kut2 = 0, Aq = — 4f2i, kd = ku, Ai = 0. 

(b) , static limit Qi/k^ = 5, for long delay kut2 = 2.0, Aq = —4^1, k^. = k^, Ai = 

(c) , fast chemical exchange (motional narrowing) Qi/k^ = 0.2, Aq = — 4^7l, Ai = 
0,kd = ku, kut2 = 0. 

Fig 3 (Color online) The 2D signal Sa signal fEqs. lfTIH) ) including TSJ and the fast spectral 
diffusion: ils/A = 0.2; ^2/^1 = 3, ka = ku = 0.002ili; Aq = -4fii, fig = A3 = Ai = 
A2 = 0; 

(a) kdt2 = 0; (b) = 1- 

Fig 4 (Color online) The 2D Sa signal (Eqs. (fTn|) ). with TSJ and slow spectral 
diffusion. ^2 /A = 5, ^2/^^! = 0.5, kd = ku = 0.002fii, Aq = -4fii, fig = A3 = 
A2 = Ai = 0; 

at times a, At2 = 0; b, At2 = 5 (i.e. kdt2 = 0.1) c, kdt2 = 1- 

Fig 5 (Color online) 2D signals for Sj generated along ki (Eq. 0), and Su generated 
along kn (Eq. Ej) 

(a) kii, At2 = 

(b) ki, At2 = 

(c) kji, At2 = 5 

(d) ki, At2 = 5 

Other parameters same as in Fig 4. 

Fig 6 (Color online) 2D signal Sa (Eqs. ^TU^ ) with spectral diffusion and TSJ. Spectral 
diffusion linewidth is varied (A2 7^ 0). Slow spectral diffusion in eg state, fast diffusion 
in overtone ef: 0.2/ A = 2, fi2/l^i = 0.6, kd = ku = 0.002fii, Aq = -4fii, A2 = 
-0.75fi2, fi3 = A3 = Ai = 0; 

(a) At2 = 0; (b) At2 = 15 (i.e. M2 = 0.1); (c) = 1- 



Fig 7 (Color online) Absorption lineshape for our model Eq.()22|) with different spectral 
diffusion width in the u and d states (^3 7^ 0). i7i = 0.5/s~^ A = OAps~^] 
Q2 = 0.33ps-\ Q3 = -0.07ps-^ kd = 0.125ps-\ K = O.lps, Aq = -2.0ps-\A2 = 0, 
A3 = Ai = 0; 



Fig 8 (Color online) The 2D signal Sa (Eq. ITTHl for the same parameters of Fig 7 at various 
time delays (a) ^2 = 0; (b) t2 = 2ps; (c) t2 = lOps. These spectra close resemble the 
experimental result o&. 
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